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The promise of innovative applications has triggered the development of many modern technolo¬ 
gies capable of exploiting quantum effects. But in addition to future applications, such quantum 
technologies have already provided us with the possibility of accessing quantum-mechanical scenarios 
that seemed unreachable just a few decades ago. With this spirit, in this work we show that modern 
optomechanical setups are mature enough to implement one of the most elusive models in the field 
of open system dynamics: degenerate parametric oscillation. The possibility of implementing it in 
nonlinear optical resonators was the main motivation for introducing such model in the eighties, 
which rapidly became a paradigm for the study of dissipative phase transitions whose correspond¬ 
ing spontaneously broken symmetry is discrete. However, it was found that the intrinsic multimode 
nature of optical cavities makes it impossible to experimentally study the model all the way through 
its phase transition. In contrast, here we show that this long-awaited model can be implemented in 
the motion of a mechanical object dispersively coupled to the light contained in a cavity, when the 
latter is properly driven with multi-chromatic laser light. We focus on membranes as the mechanical 
element, showing that the main signatures of the degenerate parametric oscillation model can be 
studied in state-of-the-art setups, thus opening the possibility of studying spontaneous symmetry 
breaking and enhanced metrology in one of the cleanest dissipative phase transitions. 

PACS numbers: 42.50.-p,42.65.Yj,42.50.Wk,03.65.Yz 


Introduction. The last decades have seen the birth 
of a plethora of new technologies working in the quantum 
regime, starting with the laser m, and including non¬ 
linear optics [HHIS], trapped ions p^HTS] and atoms [m- 
HZ], cavity quantum electrodynamics [281131] . or, more 
recently, superconducting circuits [32II36] and optome¬ 
chanical resonators [37H39] . Apart from their potential 
for quantum computation Hong and simulation [43ti5Q] , 
quantum metrology EH Eg, and quantum communica¬ 
tion [53H55] . all these technologies have allowed us to 
reach physical scenarios that were nothing but a dream 
(or a ‘gedanken’ experiment) for the founding fathers of 
quantum mechanics. 

In this work we keep deepening into the possibility 
of using new technologies to access phenomena that, 
even though predicted and theoretically analyzed since 
decades ago, have eluded observation so far, or only until 
very recently [56]. In particular, we show how modern op¬ 
tomechanical setups based on oscillating membranes [57]- 
[67] allow for the implementation of degenerate paramet- 
rie oseillation (DPO), a fundamental model in the field 
of dissipative phase transitions [68ti72] . Together with 
the laser, DPO has possibly the best-studied quantum- 
optical dissipative model, since it holds the paradigm of 
a phase transition whose associated spontaneously bro¬ 
ken symmetry is discrete ED Eg (in contrast to that of 
the laser, which is continuous ESI EH). Even though the 
main motivation for studying such a model came from 
nonlinear optics during the eighties Eg, in particular 
from the possibility of implementing it in optical para¬ 
metric oscillators, see Fig. the intrinsic multi-mode na¬ 
ture of optical cavities prevents its implementation above 


the phase transition, as we explain in the next section. 
In other words, despite the great deal of work invested 
on this model and its optical implementation, degener¬ 
ate optieal parametric oscillators (DOPOs) do not exist 
in reality. 

The situation is rather different in the microwave realm 
of electronic circuits, where one can build single-mode 
cavities in the form of simple LC circuits. Indeed, it is in 
this context where DPO has been traditionally studied in 
more detail all the way through its phase transition Est 
EH- However, an electronic circuit at room temperature 
has a very strong thermal microwave background which, 
together with other sources of technical noise, completely 
masks quantum noise and hence any possibility of analyz¬ 
ing quantum mechanical effects. This scenario was radi¬ 
cally changed with the advent of superconducting circuits 
[321136] . which are cooled down to mK temperatures, ef¬ 
fectively removing the thermal background and making 
it possible to access the quantum regime. It is in this sce¬ 
nario where, just a few months ago, quantum mechanical 
effects appearing as one crosses the phase transition of 
the DPO model have been finally observed [56] . 

Apart from being a clean system where studying fun¬ 
damental questions related to spontaneous symmetry 
breaking and ergodicity of open quantum systems [78^ 
|8T] . DPO might serve as a perfect test bed for enhanced 
metrology via dissipative phase transitions [821485] . Mo¬ 
tivated then by the interest that this model generates 
on different communities ranging from the purely theo¬ 
retical to the most applied ones, in this work we show 
that DPO can be implemented in the motion of a mem¬ 
brane dispersively coupled to the field of an optical cavity 
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FIG. 1. (a) Sketch of the optical implementation of the DPO model, (b) Phase transition of the DPO model. In the central 
panel we plot the steady-state field amplitude as a function of a, as predicted by the classical theory. The phase transition 
is revealed as non-analytic behaviour of the amplitude at threshold, and the presence of two stable solutions above threshold 
connected by the symmetry transformation a —a. The left and right panels show the steady-state Wigner function for 
a = 0.9 and 2, respectively, which for our purposes here can be interpreted as the joint probability density function W{x,p) 
providing the statistics of ‘position’ (x = a) and ‘momentum’ (p = ia^ — ia) measurements. It can be appreciated that even 
though the quantum steady state is unique and preserves the sign symmetry for any cr, it does so in two very distinct ways above 
or below threshold: with a single state centered in phase space in the first case, and with a mixture of two symmetry-breaking 
states above threshold, which deep into the phase-transition tend to the coherent states predicted by the classical limit. 


|57H63, when a multi-chromatic laser properly drives the 
latter. Starting from a first principles model, analytical 
and numerical methods allow us to identify the regimes 
where the desired model appears, as well as proving the 
feasibility of the scheme for the parameters under which 
current experiments take place. 

Degenerate parametric oscillation and its opti¬ 
cal implementation. Let us start by introducing the 
DPO model, taking its optical implementation as the 
guiding context (Fig. [^, and discussing some of the 
physics derived from it. In its minimal formulation a 
DOPO consists of an optical cavity containing a crys¬ 
tal with second order nonlinearity, and pumped by a 
laser at frequency 2uJo for which the cavity is transpar¬ 
ent (non-resonant pump configuration). The parametric 
down-conversion process occurring inside the crystal 
is able to generate photons at the subharmonic frequency 
ujo, assumed resonant. The model can then be formu¬ 
lated as the following master equation for the state p of 
the intracavity field [86l [87] : 


dp 

dt 


^DP0 5 P 


l^VAp]+lVa[p\, ( 1 ) 


with 


ffopo = woaU + (2) 


where d is the annihilation operator of cavity photons, 
and we use the notation T>j[p] = 2JpJ^ — Jp — pJ^ J. 
The last term describes the loss of cavity photons through 
the partially transmitting mirror (with corresponding 


damping rate 7 , proportional to the mirror transmit¬ 
tance). The second term describes the loss of photon 
pairs which, after being down-converted into a pump pho¬ 
ton, leave the cavity to never come back (at rate 
where g is proportional to the crystal’s nonlinear suscep¬ 
tibility). Finally, the Hamiltonian term describes the co¬ 
herent exchange of photon pairs with the pumping field 
via down-conversion (at rate 7 cr/ 2 , where a is propor¬ 
tional to the amplitude of the laser), as well as the free 
evolution of the cavity mode. 

The first thing to note is that this master equation is 
invariant under the parity transformation IJ = (— 1 )"^^, 
which performs the operation WdU = —d. On the other 
hand, defining a = lim^^oo it is well known 

[Il|88l|89| that the classical limit of this equation pre¬ 
dicts an off (or below-threshold) stationary state o; = 0 
for cr < 1, and an on (or above-threshold) phase-bistable 
state a = ±y/ 2 (cr — l)lg for cr > 1 |89]. Hence, at cr = 1 
(threshold) the classical theory predicts a phase transi¬ 
tion, accompanied by spontaneous symmetry breaking of 
the discrete phase above threshold, since the system has 
to choose between two possible steady states which indi¬ 
vidually do not preserve the symmetry. 

In contrast to the classical state, the quantum steady- 
state solution p = limt^oo of Eq. Q 

is unique for any cr [zuiziEn]- The symmetry of the mas¬ 
ter equation, together with the uniqueness of the steady 
state, forces the latter to be invariant under the trans¬ 
formation as well, UpU^ = p, what in turn implies that 
a = tr{pd} = 0 Vcr. This could lead to the conclusion 
that quantum theory completely spoils the phase transi- 
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FIG. 2. (a) Sketch of the optomechanical setup with which we propose to implement the DPO model, (b) Asymptotic phonon 
number as a function of the pump parameter, for a set of parameters in correspondence with current experiments (see the text 
for details). The main plot corresponds to the predictions in classical limit, and the inset to the predictions obtained from the 
asymptotic state generated by the master equation. The blue circles correspond to the optomechanical model, while the red 
curves are evaluated from its effective DPO model. Note the good agreement between both models. 


tion and its associated spontaneous symmetry breaking. 
However, the situation is a bit more subtle: as shown 
in Fig. through the Wigner function [ 8 ^ fOTHOT] . be¬ 
low threshold the quantum state is a squeezed state cen¬ 
tered at the origin of phase space, while above threshold 
it develops two lobes centered (approximately) around 
the classical bi-stable solutions. This shows that in the 
quantum domain the phase transition predicted at the 
classical level has the significance of a crossover between 
phases which preserve the symmetry in two physically 
distinct ways. 

Unfortunately, in real experiments degenerate down- 
conversion has to compete with non-degenerate channels 
in which the photon pairs are generated in two different 
modes with frequencies uji and UJ 2 such that uji-\-uj 2 = 2 co’o 
{energy eonservation)^ and it is possible to show that 
phase-matching in the nonlinear crystal {momentum eon- 
servation) always gives preference to one of such pro¬ 
cesses above threshold [QSfllQQj . Therefore, these devices 
cannot be used to study experimentally the DPO model 
all the way through its phase transition. 

Optomechanical implementation of degenerate 
parametric oscillation. In contrast to the optical case, 
we show in this work that optomechanical resonators in 
which a mechanical degree of freedom is dispersively cou¬ 
pled to the cavity field allow for the implementation of 
DPO all the way through its phase transition. Our pro¬ 
posal follows closely current experimental setups based 
on dielectric membranes embedded in optical cavities 
[52H67]. In such setups, the type of coupling arising be¬ 
tween the membrane’s motion and a given driven cavity 
mode depends on the position of the former with respect 
to the standing wave defined in the cavity by the latter 
[571 [sa EH Eg. In particular, denoting by x the dis¬ 
placement of the membrane with respect to its equilib¬ 
rium position (normalized to its zero-point fluctuations 


the frequency shift felt by the optical mode is pro¬ 
portional to when the membrane is located in a node 
or an antinode of the mode’s standing wave, while it is 
proportional to x when it is half-way between them. In 
most optomechanical systems the linear coupling domi¬ 
nates, and many exciting phenomena have been already 
proven by exploiting it, including mechanical cooling 
jsziEoiEaEHiiniHins], optical squeezing [IIQHII31, and 
induced transparency [66l I114flll7j . On the other hand, 
the quadratic coupling has already promised very inter¬ 
esting applications such as quantum nondemolition mea¬ 
surements of the phonon number ISZIISH], and in our case 
it will provide the leading mechanism to achieve degen¬ 
erate parametric oscillation. 

Our proposal is sketched in Fig. lit- We consider 
two optical modes with frequencies uj\ and cjq, linearly 
and quadratically coupled to fundamental mechanical 
mode of the membrane (with frequency U), respectively; 
consequently, we will refer to them as the linear and 
quadratie modes. This can be achieved with a single 
optical cavity by selecting two different resonances with 
the appropriate relative separation between nodes and 
antinodes. In the absence of optomechanical coupling, 
the membrane is at some equilibrium temperature T 
with its thermal environment, which drives it at some 
rate to a thermal state with mean phonon number 
^th = ksT/hfl. However, we assume the linear mode 
to be driven by a monochromatic laser at frequency cjl 
tuned to cool down the membrane to an effective phonon 
number fieff ~ nth/C\ at rate 7 eff = Ci7m, where C\ is 
the cooperativity of the linear optomechanical coupling 
pUTlfTOH] : current experiments reach cooperativities on 
the order of 10000 , which together with cryogenic tem¬ 
peratures allow to cool down the mechanical motion close 
to its ground state m ^eff = O 5 what we assume in the 
following. On the other hand, the quadratic mode is 
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driven with a laser containing a tone at frequency cjq, 
plus a sideband at frequency ujq — ftq. We will use the 
first tone to create the two-phonon losses needed in the 
DPO model Q, while the combined action of the two 
tones will provide the coherent exchange of phonon pairs. 

We model the system by a master equation governing 
the evolution of its state p, which in a frame rotating at 
the laser frequency ujq takes the form 

P]+7q^aq[p]+7eff^6[p]5 (3) 

with Hamiltonian terms 

= Qpb, (4a) 

Hq = Aqd^dq + i[fq(t)a^ - (t)aq], (4b) 

-^qm ^q^ ^q^q: 

where the bichromatic driving amplitude of the quadratic 
mode can be written as Sq{t) = fq + ^ with <p 

some relative phase between the two tones. Aq = ujQ—ujq 
is the laser detuning, h is the mechanical annihilation op¬ 
erator, from which the mechanical displacement is writ¬ 
ten as X = 6 - 1 - 6 ^. dq is the quadratic mode’s an¬ 
nihilation operator, with corresponding cavity damping 
rate 7 q (proportional to the mirror transmissivity). The 
driving amplitudes fo,i can be written in terms of the 
power Pop of the laser at the corresponding frequency as 
5o,l ~ \/27qPo,l/dcJq [HH]- 

We can understand the conditions under which this 
bichromatically-driven optomechanical model is mapped 
to the DPO model by adiabatically eliminating the 
quadratic mode. We provide the details of the deriva¬ 
tion in [89], and here we just want to point out some 
physically relevant steps. We follow the usual projector- 
superoperator technique |118I4121] in which the quadratic 
mode is assumed to be in some reference state and follow 
some reference dynamics. For our current purposes, it is 
enough to assume that it does not feel any mechanical 
backaction, so that its dynamics is described by the mas¬ 
ter equation above with ^q = 0. This means that: (i) we 
can take a coherent state with amplitude 

aq{t) = (5) 

as its reference state, where we have made a concrete 
choice of the second tone’s phase 0 that simplifies the 
expression [89], and defined fij = Ej/["fq + (Aq + jif)q)^], 
which are interpreted as the number of photons intro¬ 
duced by the corresponding laser in the cavity; and (ii) 
all the correlation functions of its quantum fluctuations 
6aq = aq — aq will decay in time at rate jq. 

Keeping in mind the traditional picture of sideband 
cooling |lQ7l41Q9j . it is intuitive to understand how the 
DPO model arises from Hqm upon adiabatic elimination 
of the optical mode. First, it will turn out to be conve¬ 
nient to work in the weak sideband regime ni ^ no [HS] . 


The coherent part of the optical field generates then an 
effective mechanical Hamiltonian 

H,s = ^b^b-gMt)\^3:^ ( 6 ) 

« flesPb + igq\/noni(e‘^‘'*S^ — H.c.), 

where Oeff = ft — 2 ^qno, and any other term can be 
neglected as long as we work within the rotating-wave 
approximation ^ gqfio jHl]. Provided that the 

second tone is chosen as ftq = 2 f 2 eff, this effective 
Hamiltonian provides precisely ^dpo as introduced in 
(H). On the other hand, the elimination of the opti¬ 
cal fluctuations generates both dissipative and Hamil¬ 
tonian mechanical terms. Within the weak sideband 
and rotating-wave approximations, all the Hamiltonian 
terms can be neglected, while only the dissipators P 52 , 
1 ) 512 , and survive, corresponding, respectively, to 

two-phonon cooling, heating, and dephasing. In the 
weak sideband regime, the rate of each process is static 
and solely controlled by the fundamental tone. Hence, 
setting its detuning to the red two-phonon resonance, 
Aq = — 2 f 2 eff, while working in the resolved sideband 
regime, 4^4^^^ ^ 7 q, the heating and dephasing terms are 
highly suppressed [89] just as in standard cooling |1Q7F 
ung, leaving us with a two-photon cooling dissipator at 
rate ymC'q, where we have introduced the cooperativ- 
ity Cq = gqfio/iqlm- H IS important to note that the 
Markov approximation (central to this method), requires 
a decay of the optical correlators much faster than the ef¬ 
fective mechanical dynamics induced by the optical fields 

inn], that is, 7q » max{7eff,7mC'q,flqVh^}. 

This shows that the dynamics of the mechanical mode 
should follow an effective DPO master equation of the 
form with = fieff, 7 = 7eff, 9^ = 4Cq/Ci, and 

= ( 7 \/h.i 7 q/ 7 eff. Remarkably, we obtain a set of pa- 
rameters that can be optically tuned by means of the 
two laser powers Pq and Pi and that, as we will show 
below, allow to explore physical regimes not available in 
all-optical implementations. In the following we discuss 
whether the parameters which can be reached in current 
experimental setups, together with the bounds that the 
above-mentioned conditions impose on these, are com¬ 
patible with having a reasonable range for these effective 
DPO parameters. 

Implementability in current setups. We take the 
experiments of m as a reference, in which Q = 4.4MHz, 
7 m = 0.8Hz, cjq = 1770THZ, 7 q = 1.3MHz, Q = 10000, 
and gq = 10“^Hz (this last parameter taken from [65]). 
In order to stay safely within the rotating-wave approxi¬ 
mation we must impose a bound on the intracavity pho¬ 
ton number given by nq ^ ft/^^gq ~ 10 ^^, what also 
certifies that we are within the resolved sideband regime, 
^eff/ 7 q ~ 3.5. Taking nq = 3.3 x 10^ (corresponding to a 
power of Pq ~ 20mW, fairly reasonable in such setups), 
we then obtain a quadratic cooperativity Cq ~ 3.3 x 10“^, 
leading to an effective two-phonon loss g ~ 10 “^, on the 
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order of the one obtained in optical implementations [88] . 
On the other hand, let’s take ni = 3 x 10^ ^ no (cor¬ 
responding to Pi ^ 35/iW) as an upper bound for the 
sideband photon number; then by varying the sideband 
power from zero to this value, the effective a parameter 
can be varied all the way through the phase transition 
and up to a = 2.5, showing how current optomechanical 
setups should be able to reach regions of the DPO model 
beyond what’s possible in optical implementations. Note 
finally that the Markov approximation is very well satis¬ 
fied since 7q/7eff ~ 150. 

Numerical simulations. In order to certify the 
predictions offered above with the effective mechanical 
model under optical adiabatic elimination, we have per¬ 
formed numerical simulations of the full optomechani¬ 
cal problem for the realistic parameters of the previous 
section. We proceed in two ways. First, by projecting 
the master equation (|^ in a truncated Fock basis for 
the optical and mechanical modes, allowing us to di¬ 
rectly simulate the evolution starting from any initial 
state (89] [91]. Since the master equation is manifestly 
time-dependent, so will be the asymptotic state; in par¬ 
ticular, we are interested in the asymptotic phonon num¬ 
ber which oscillates at frequency flq and 

can be approximately written as n + (5nsin(flqt), where 
typically n ^ Sn. In the inset of Fig. we show the 
static phonon background n as a function of the effective 
a (which can be tuned through the sideband power Pi 
keeping the rest of parameters fixed), together with the 
phonon number predicted from the effective DPO mas¬ 
ter equation §■ For large phonon numbers this type 
of “brute force” simulation becomes unfeasible, and does 
not allow us to get above a ~ 0.95. Fortunately, in order 
to prove that the DPO phase transition is present in the 
full model, it is enough to consider the classical limit of 
the model, which is the second type of simulation that 
we have performed. In this limit the optical and me¬ 
chanical modes are described by complex amplitudes a 
and respectively, which, defining the mechanical po¬ 
sition X = 2Re{;d} and momentum p = 2Im{/d}, evolve 
according to [89] 

X = Up, (7a) 

p = -27effp -{n- 4pq|(a|^)x, (7b) 

= ^q(t) (7q i^q )Q^- 

These are a set of coupled nonlinear equations which can 
be efficiently simulated in virtually all parameter space. 
They predict an asymptotic phonon number given by 
limt^oo[^^(^) +P^(^)]/4, the static part of which we plot 
as a function of a in Fig. We also show in the figure 
the predictions of the DPO model in the classical limit, 
which are analytical and given hy n = 2(a — l)/g^. We 
see that both the quantum and classical simulations find 
very good agreement with the DPO model, in particu¬ 
lar the classical limit, which shows the phase transition 


exactly as expected. 

Conclusions. In summary, we have shown that the 
elusive degenerate parametric oscillator model can be re¬ 
alistically implemented in current optomechanical setups. 
Apart from providing the possibility of studying experi¬ 
mentally many interesting theoretical predictions put for¬ 
ward during the last three decades, the implementation 
of this simple (but paradigmatic) dissipative model in 
modern quantum technologies opens the way to analyz¬ 
ing open questions related to ergodicity and spontaneous 
symmetry breaking, as well as enhanced metrology with 
dissipative phase transitions. 
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Supplemental material 


This supplemental material is divided in three sections. In the first one, we explain how we obtain the asymptotic 
state of the degenerate parametric oscillator (DPO) model, both in the quantum and classical regimes. In the second 
section we make a detailed derivation of the effective mechanical master equation under adiabatic elimination of the 
optical mode, stressing the conditions and approximations under which it is expected to be equivalent to the DPO 
model. In the final section we explain how we have found the asymptotic state of the complete optomechanical model, 
again both in the quantum and classical regimes. 


ASYMPTOTIC STATE OF THE DEGENERATE PARAMETRIC OSCILLATOR 

Moving to a simpler picture. As introduced in the main text, the master equation corresponding to DPO 
can be formulated for a single bosonic mode (corresponding to the down-converted intracavity mode in the optical 
implementation) and it reads 


dp 

dt 
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Va2[p] ^ -fVa[p], 


( 8 ) 


with 


^DPo(t) = ujoa^d + (9) 

where we remind that d is the bosonic annihilation operator, and we use the notation T>j[p] = 2JpJ^ — Jp — pJ^ J. 

The parameters are all defined in the main text. 

The asymptotic state of the DPO is better analyzed in a picture rotating at frequency cjq, where the state becomes 
time-independent. In particular, defining the transformation operator Uc = exp(icc;o^d’*'d), in the new picture the state 
of the system p = UcpUl obeys the master equation 

2 

+ [p] + 7 ^a [p] , (10) 

where the transformed Hamiltonian reads 

^DPO = ^c^DPO^c “ ^od‘*'d = i7cr(d‘*‘^ - d?)/2. (11) 

Hence, we see how in this picture the master equation becomes time-independent, leading to a stationary asymptotic 
state of the system. 

Quantum asymptotic state. Indeed, the unique steady state of this master equation is known analytically, in 
particular in the form of a positive P distribution m, from which in principle the elements of the density operator 
can be reconstructed in any basis |122j . However, such a reconstruction is computationally very demanding, and for 
our purposes it is simpler to evaluate the steady state numerically. Concretely, in the main text we have shown the 
Wigner function associated to the steady state for ^ = 0.1 and two different values of the pump parameter, a = 0.9 
and <7 = 2 , above and below the phase transition, respectively. Let us now spend some time explaining how we have 
computed this steady states and their corresponding Wigner functions. 

As explained in detail in m, we perform the numerics by going to superspace, where the elements of the den¬ 
sity matrix in the Fock basis are gathered in a vector p, and the master equation becomes then a linear system 
dp/dt = LdpoP, where Ldpo is then a representation of the Liouville superoperator which induces the DPO dynam¬ 
ics, >Cdpo['] = ~i[^DPO:'] +7^a['] + 7^^^a2[']/4. Notc that the Fock basis is infinite-dimensional, and hence one 
has to introduce a truncation {|n)}n=o,i,...,Ar in order to work in the computer. In our case, we follow the criterion 
of truncating to values of N for which the observables we are interested in (e.g., the photon number) converge up to 
a three-digit precision or more. The steady state corresponds then to the eigenvector with zero eigenvalue of Ldpo, 
which provides us with the Fock basis components of the steady-state operator, {pmn}m,n=o,i,...,Ar- 

Once we have the density operator in the Fock basis, the Wigner function can be evaluated as follows. First, given 
a harmonic oscillator with position and momentum quadratures, x = d) Pd and p = i{d^ — a), respectively, recall that 
the Wigner function W{x,p) can be seen as a joint probability density function for measurements of such observables, 
in the sense that the marginal P{x) = J^dpW{x,p) provides the probability density function predicting the statistics 
of position measurements (and similarly for the momentum). Defining the polar coordinates (r, (p) in phase space by 


ap 

d.t 


^DPO, P 
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(x^p) = r{cosp, sin(/9), the Wigner function of the steady state p can be found from its components in the Fock basis 
as 


N 


W{r,ip)= y] PmnWmn{r,^), 


mn=0 


where we have defined the Wigner function of the operator |m)(n|, given by 


Wmn{r,(p) = 


i-iy 


' ^m-nJ^m-n/2^ 


( 12 ) 


(13) 


with L^{x) the modified Laguerre polynomials and where we have assumed m> n (note that Wnm = ^mn)- 

Classical limit and steady state. Phase transitions in dissipative systems are usually revealed in the classical 
limit of the corresponding quantum models. Let us explain how such limit can be obtained from the master equation 
describing the system quantum mechanically. The idea is indeed quite simple in the case of bosonic systems: the 
classical limit consists in assuming that the state of all bosonic modes is coherent, with an amplitude that will play 
the role of the classical variable. 

In particular, in the case of the single-bosonic mode considered in the DPO, this means that we assume its state to 
be a coherent state \a{t)) at all times, such that the expectation value of any normally ordered observable factorizes 
as {t)d^(t)) = ^ where for convenience we are defining expectation values with respect to the state in 

the picture rotating at the laser frequency, that is {t)d}{t)) = tr{d’*'^d^p(t)}. This allows us to find an evolution 
equation for a{t) from the master equation (10) as follows: given any operator A, the master equation allows us to 
write the evolution of its expectation value as 


d{A) 

dt 


= trM 


dt 


= -i([i, Hnopo]) + AYa^) + «'])) + 7 (([at a])), ( 14 ) 


which applied to the annihilation operator, and using the coherent-state ansatz provides us with the classical equation 
of the DPO 


-1 • * 9 \ \2 

7 a = aa —a — a. 


(15) 


This is indeed the equation that would have been obtained by using classical electromagnetic theory on the DOPO, 
where a would be interpreted as the normalized amplitude of the optical field. As explained in the text, this 
equation has two types of asymptotic, stationary {a = 0) solutions: the trivial one d = 0, and a nontrivial one 
a = ±^2((7 — 1)/^, which exists only for a > 1 and has sign-indeterminacy owed to the symmetry a ^ —a oi Eq. 

We use the bar to denote “stationary state” In order for these solutions to be physical, they need to be stable 
against perturbations; their stability can be analyzed by studying the evolution of small perturbations around them. 


that is, by writing a{t) = a-\-6a{t)^ and linearizing Eq. (15) with respect to 6a. Defining the vector 6cx = col((5(a, (5a*), 
one obtains the linear system dScx/dt = AiScx, where the linear stability matrix reads 


M = 


-1 - 2g 




a-g^ay2 


a — g‘^a'^^12 —l — 2g‘^\a\‘^ 


(16) 


Hence, the stability of a given stationary solution a is determined by the eigenvalues of this matrix: when they all 
have negative real part, it will be stable, while if some of them have positive real part, perturbations will tend to 
grow, showing that the solution is unstable. In the case of the trivial solution, the eigenvalues are A± = — 1 ± a, 
and hence, it is unstable for a > 1. On the other hand, the eigenvalues associated to the nontrivial solution read 
A± = —2(7 + 1 ± 1, which are always negative for a > 1, and hence this solution is stable. 

Therefore, we see that at the classical level the phase transition is revealed by a non-analytic change in the stationary 
solution at threshold (7 = 1. 


ADIABATIC ELIMINATION OF THE OPTICAL MODE 


Moving to a simpler picture. Our starting point is the master equation of the optomechanical system as we 
introduced it in the main text: 

= —i[^m + ^q(^) + ^qm, p] + lq^a^[p] + 7eff ^6 [p], 


(17) 
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with Hamiltonian terms 


Hn, = nh^b, 

(18a) 

Hq = AqaJjfflq + i[£’q(i)a)j - £qit)aq], 

(18b) 

Hqm ~ ^7q^7q^ ’ 

(18c) 


and where the bichromatic driving amplitude of the quadratic mode can be written as fq(t) = ^ with 

(j) some relative phase between the two tones that will be chosen shortly. All the symbols have the meaning introduced 
in the main text. 

In order to perform the adiabatic elimination of the optical mode, it is convenient to move to a picture where the 
large coherent background that the driving fields create in the optical mode is already taken into account. This is 
accomplished by using a displacement D[ac^{t)] = exp[(aq(t)aj — (a*(t)dq] as the transformation operator, where the 
amplitude (aq(t) is chosen to obey the evolution equation 


dq — ^q(^) (Tq 


(19) 


with solution 


a. 


,(i) = + 


fo 


7q - iAq 


I _ g-(7q-iAq)t 




7q i(^q “I" ^q) 


g _ g (7q iAq)t 


( 20 ) 


Note that choosing 0 = —7r/2 + arctan(Aq/ 7 q) — arctan[(Aq + flq)/ 7 q], we obtain the asymptotic displacement 


lim aq{t) = 

t>7q ^ 


( 21 ) 


with no = ^o/[7q + ^q] ~ ^i/[7q + (^q + ^q)^]^ which is the amplitude that we introduced in the main text. 

In the following we assume to be working in this asymptotic regime t ^ 7q"^, even if we don’t write the limit explicitly 
to shorten the expressions. In this new picture, the transformed state p = &[ac^{t)]pD[ac^{t)] evolves then according 
to 


^ p] + 'JqDa^lp] + '^eS'Dblp] 


( 22 ) 


where the transformed Hamiltonian 


H{t) = D^[aq{t)][Hm + Hq{t) + Hqm]D[aq{t)] + i (dqfflq - dqfl^) , (23) 

can be written as the sum of three terms, H = + H^(t) + i7qm(0, with 

= Aqdqhq, (24a) 

H^{t)=nh^h-gq\aq{t)\^x\ (24b) 

ffqm(i) = -gq[al{t)aq + aq{t)al]x^ - gqal&qx'^. (24c) 


The interest of moving to this picture is that now the driving has been moved to the coupling and the mechanical 
Hamiltonians, what will allow us to easily understand the physics behind the system. Moreover, in the absence of 

— i[Aq( 2 q( 2 n 


qL J ~ H^q^q^q^ J ~ 
which in the original 


optomechanical coupling the dynamics of the optical mode is generated by the Liouvillian £, 

7 qPaq[-], including only detuning and dissipation, which drives it to a vacuum state at rate 7 q 
picture corresponds to a coherent state with amplitude a^{t )—. 

Derivation of the effective mechanical master equation. Let us rewrite the master equation (22) as 


f^=C<^)[p]+Cq[p]+C^l[p], 


(25) 


where £q is defined above, while £m^[-] = —i[^ni(^)r] +7eff^6['] and >Cqm['] = —i[^qm (^)5 *]• Adiabatic elimination 
proceeds by choosing some reference state and dynamics for the optical mode, which is assumed to remain unperturbed 
by the mechanical mode, and hence the accuracy of the elimination depends crucially on the choice of a proper 
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reference. For our purposes, it is enough to take the dynamics generated by >Cq as the optical reference, and hence its 
steady state pq = | 0 )q( 0 | (vacuum) as the reference state. Let us then define the projector superoperator 


V[-] = Pq®trq{-}, 


(26) 


and its complement Q = 1 — V. These superoperators satisfy the useful relations = CmV^] (obvious since 

Cm acts on the mechanics only) and VCq[-] = 0 = Tq'P[-] (where the second equality is again obvious, while the first 
one comes from £q being traceless by conservation of probability). 

The next step consists on projecting the master equation onto the corresponding subspaces defined by these super¬ 
operators. Defining the projected components of the density operator u(t) = V[p{t)] and w{t) = Q[p(t)], and using 
the properties of the projectors, it is straightforward to get the coupled linear system 


^ = [u]+VC^*i[w], 

= (-^m + A + ^4m) + 24m N- 

The second equation can be formally integrated, leading to 



(27a) 

(27b) 


(28) 


where T is the time-ordering superoperator, and we have not written the term which depends on the initial value w( 0 ) 
since in this concrete dissipative scenario any information related to it has to be completely washed out asymptotically. 
Next, we introduce this formal solution in the first equation, and perform a Born approximation in which we neglect 
terms beyond quadratic order in the interaction >Cqm, what allows us to write 


du 

dt 


= £W[u]+p44u] + 


Jo 




(29) 


where in addition we have made the variable change t' = t — r in the integral, used the fact that >Cm and £q commute, 
and defined the mechanical time evolution superoperator 


^ 7 - |e/*-x I ^ (30) 

After performing the partial trace over the optical mode, this equation provides an effective master equation for 
the reduced mechanical state pm = trq{p}. In order to simplify further such equation, we need to write the explicit 
form of the interaction >Cqm, what we do as 


3 

i=i 

with 

B = (Uq , 4, ^ = flq [“q ’ “q a] • 


(31) 


(32) 


Note the null asymptotic expectation value of the optical operators, trq{ 5 jPq} = 0 Vj, meaning that the second term 
of Eq. (29) does not contribute since VC^m — 0- Let us then define the optical correlation functions 


tVq{Bie^<^^[pqBj]} = lim {Bj{t)Bi{t + T))q = Kjiir), (33a) 

t^OO 

trq{Bie^-^'^[Bjpq]} = lim {Bi{t + T)Bj{t))q = Hjiij), (33b) 


where the expectation value refers to the picture rotating at the laser frequency and we have used the quantum 
regression theorem |12Qj . Then, the effective mechanical master equation can be written as 


dt 


^ -^m [Pm] + T Gj{t)Gl{t) f dTKji{T){x'^U^~'"'*'>[pm{t - t)x'^] - [pm{t “ t)x'^]x'^} (34) 

jl=l “'o 

N .t 

W Gj{t)Gl{t) / dTHji{T){U^~'^'*^[x^Pm{t - t)]x^ - Pm{t “ "t)]}. 

ii=i -^0 
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The correlation functions Kji{T) and Hji{T) can be evaluated in many ways. Instead of using the dynamics induced 
by >Cq in the Schrodinger picture, a particularly simple way of evaluating them is by using the equivalent quantum 
Langevin equations of the optical operators, which in the simple case of having dissipation and detuning only, consist 
of a single closed equation for the annihilation operator |119j : 


ddq 

dt 


(Tq i^q)^q T ^ 


where the only non-zero input-operator correlators up to fourth order are 

(ain(^i)aL(^2)) = S{ti - h), (ain(il)4n(^2)ain(^3)a|„(^4)) = - t2)5{t3 - q), 

{ain{ti)ainit 2 )al^{t 3 )al^{t 4 ,)) = S{ti - t 3 )S{t 2 - U) + 5 {ti - U) 5 {t 2 - ^ 3 )- 

The asymptotic (t ^ 7q"^) solution of this equation reads 

«q(i) = / dt'ain{t'), 


(35) 

(36) 


which together with the correlators of the input operator allows us to write 

Kji{T) = 1(5,2, and if,-,(r) = 2(5,1, 


(37) 


(38) 


which are functions decaying at rate 7 q. 

Hence, of the dynamics induced by the time evolution superoperator Z7m we see that only the processes 
happening at a rate faster than or similar to 7 q play a role in the integral terms of the effective mechanical master 
equation (34). This brings us to the final major approximation, known as the Markov approximation: we assume 
that of all the processes contributing to the mechanical dynamics, the only term acting appreciably on the time- 
scale of the optical decay is the simple oscillation induced by the term [ft — 2gc^{no + ni)]b^b of Hm (we will justify 
this approximation self-consistently at the end of the derivation). Within this Markov approximation, we can then 
approximate 


- t)x^] ~ « Pm{t)x{Tf 

where fleff = ^ — 2^q(no + tii) and x{t) = and similarly 

- t)] x{T)^Pm{t), 

leading to the effective mechanical master equation 

+ {[r(t; + r(t; + r(t; 0){2b^b + l)]pn 


dpva 

dt 


(39) 

(40) 

(41) 


—X 


r(i; + r(i; + r(i; 0 ){ 2 b^b+ 1) 


Pm + H 


•C.}, 


where, after performing the time integrals (in the asymptotic limit), we have obtained asymptotic time-dependent 
rates 


r(t;w) = 


TmCq 


1 - i(Aq + 2w)/7c 


(^1 + iyni/noe**^"*) + 


7m Gq 


— i(Aq + Q,q 


■ 2w)/7c 


(ni/no - iyni/noe , (42) 


with the cooperativity defined as Cq = 5'qno/7q7m- 

In order to see that this master equation has all the ingredients that we need plus many more, so it is just a matter 
of finding the regime in which the latter do not contribute, let us rewrite it. By defining the real and imaginary parts 
of the rates, rR(t;ci;) = Re{r(t;ci;)} and ri(t;co’) = Im{r(t; u;)}, and defining the phonon-number operator n = b'h, 
we can write 


dpm 

dt 


Pjji +7eff^6[Pm] +rR(^; ^eff )^62 [Pm] +rR(^; ~^eff )^6t2 [Pm] +4rK(t; 0)Pn[Pm] +>^NRw[^m]5 (43) 


where we have defined the effective Hamiltonian ^eff(^) = ^dpo(^) + ^±dpo (^)5 containing terms that we will need 
for the DPO model 


Hi)Po{t) = flesfi + igq^/nonl{e 


(44) 
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plus some that we don’t want to contribute 

H±DPoit) = ['^gq^/f^onl sin{Qqt) -Ti{t; Qes) + STiit; -Oefi) + 4ri(t; 0)]n (45) 

+ fieff) + ri(i; -f2eff) + 4ri(t; 0)]n^ - [flq(no + ni + + H.c.], 

and we have collected into £nrw the terms which in the absence of the sideband are expected not to contribute within 
the rotating wave approximation, which read 

^RwIPm] = [r(t; fieff) + r*(i; -n,s)ppj‘^ - T{t-, - T* {t; -n,s)pj>‘^ (46) 

+ [r(i;^^eff) + r*(t; 0)]&^Pin(2n+ 1) - r(i; fieff)(2n + l)Ppyn - T*{t\Q)pya{2h + l)P 
+ [r{t; 0) + -l^eff)](2n + l)/5m6" - r(t; 0)b‘^{2h + l)p^ - r*{t; -n^s)pmb\2h + 1) + H.c.. 


Degenerate parametric oscillation regime. Let’s now discuss the conditions under which the effective me¬ 
chanical master equation above will correspond to the master equation of the DPO, Eq. (§. 

Looking at the term L^dpo {t) of the effective Hamiltonian, we see that the sideband flq should be chosen to match 
twice the effective mechanical frequency, that is, 

On the other hand, we would like the effective two-phonon cooling P 52 to dominate over any other irreversible 
process, in particular over the two-phonon heating I> 5 t 2 and dephasing and to do so with a time-independent 
rate. Looking at (42), the latter can be naturally accomplished by driving the fundamental tone much stronger than 
the sideband, that is, no ^ fii- The static part of the rates (42) then suggests that cooling will be enhanced by 
choosing a detuning of the fundamental driving tone matching the red two-phonon sideband, Aq = — 2 Peff; indeed, 
this choice provides the following real, static part of the rates 


Tcooling r];ieating 


TmC'q 


1 + ’ 


and r, 


dephasing 


IrnCq 




(47) 


showing in addition that we need to work in the resolved sideband regime ^ 7 ^ in order for heating and 

dephasing to be suppressed; in particular, we will define the parameter r = (1 + 1 , which allows us to 

approximate the rates by 


r(t; f2eff) = 7mC'q (^1 + i\/rni/no + + i/r(ni/no)e“^‘*^'*"*^ , (48a) 

T{t; -fieff) = -iVrimCq (^1 + - 2i^ni/noe~^'^^^^^^ /2, (48b) 

T{t; 0) = 7 mC'q (fii/ho - , (48c) 


expressions that together with working in the weak sideband regime, fiilno <C 1 , suggest that the only relevant rate 
is the real static part of r(t; Peff) which provides the two-phonon cooling rate as desired. 

The next constrain on the parameters comes from the fact that there are many counter-rotating terms that we 
want not to contribute within the rotating-wave approximation. Among these, inspection of and iL^DPO shows 

that the largest of such rates are ymC'q and gq{no + ~ ^qf^o, but note that ^qno/TmC'q = 7 q/^q which is typically 

much larger than 1 (optomechanical systems work far from the single-photon strong coupling regime, specially when 
referring to the quadratic coupling), and hence ^qU-o is the largest of these two rates. Thus, the rotating wave- 
approximation requires ^qU-o ^ Provided that this approximation holds, we can neglect all the counter-rotating 
terms in ^xdpo and approximating them by 


^±DPO - -v^ 7 mC'q(ll + 9h)n/2, 


(49) 


and 

^NRw[/5m] = y^i 7 ^ 7 mC'qe^‘^"“*[iS^/ 3 m( 2 n + 1) - i(2n + l)b‘^pm + Vrp^{2h + 1)V^ (50) 

- y/rP{2h + l)pra + Vrpmb‘^{2n + 1)] + H.c.. 


This expression clearly shows that T^rw negligible when compared with the two-phonon cooling term 7 i„( 7 qDf, 2 . 
On the other hand, ^±dpo provides a negligible effective mechanical frequency shift, but also a Kerr term which is 
expected to be negligible only as long as (n) ^ f^eff/d.SyTymC'q, which puts a bound on the number of phonons. 
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Nevertheless, for the parameters corresponding to a realistic implementation used in the main text, we find this bound 
to be ^ 10^^, while the classical limit of the DPO tells us that the phonon number expected at cr = 2.5 is on the order 
of 10^^, three orders of magnitude below the limit in which the Kerr term can start playing a role. 

Provided that all these considerations are taken into account, we then expect the effective mechanical master 
equation to be very well approximated by 


with effective Hamiltonian 


dpiYi 

dt 


(^): Pn 


'Teff^6[Pm] [Pm]: 


(51) 


ffeff(i) ~ + ^9qVn(^{e 2i«ettt5t2 _ g2inefft^2^^ ^52) 

which is exactly the DPO model, Eq. §. 

Note that there is a final constrain on the parameters coming from the Markov approximation that we performed 
in the adiabatic elimination: since we assumed that within the decay of the optical correlators at rate 7 q the only 
relevant mechanical process is the simple oscillation at frequency we need the rates 7eff, 7 mC'q 5 and ^qv^noPi, 
to be smaller than 7 q. For the parameters considered in the main text, 7eff is the largest rate of the three (for an ni 
corresponding to a = 2.5 in the DPO model or smaller), and it satisfies 7eff/7q ~ 6 x 10“^, so we are safely within 
the Markov regime. 


ASYMPTOTIC STATE OF THE FULL MODEL 


Let us in this final section of the supplemental material explain how we have found the asymptotic state of the 
optomechanical model numerically. As with the stationary state of the DPO model, we have performed simulations 
of the full master equation ( pT| ), as well as simulations in the classical limit. 

In the first case, our starting point has been the optomechanical master equation in the displaced picture, Eq. 


(22). For numerical purposes, it is important to work in this displaced picture because the state of the optical mode 


should stay close to vacuum, while in the original picture it is a highly populated coherent state which does not 
allow for a reasonable truncation of the optical Fock space. As before, we set the truncation of the mechanical and 
optical Fock bases in such a way that the mean phonon number finds convergence up to the third significant digit, 
what typically does not require more than one or two photons in this displaced picture. The simulation proceeds 
again as explained in detail in [91], that is, by moving to superspace where the master equation (22) is turned into 


a linear system dp{t)ldt = L(t)p(t). Note that now the linear problem is manifestly time-dependent, and therefore, 
there will be no steady state. In particular, L(t) is 27r/Uq-periodic, and this periodicity is reflected into a time- 
dependent asymptotic state p{t) = which we find by solving the linear system numerically starting 

from different initial conditions p(0), that is, different initial states p(0). In all the simulations we have checked 
that the asymptotic state is independent of the chosen initial state (e.g., vacuum or the steady state of the DPO 
for the mechanics). As explained in the text, the observable we have focused on is the asymptotic phonon number 
= tT{b^bp{t)}, whose time evolution can be approximated by a function of the type n + (5nsin(Uqt), 

with Sn n. 

The superspace simulation of the master equation becomes quite heavy as the mechanical state gets populated, 
what has prevented us from performing simulations for values of the sideband power where the system is expected to 
be above the DPO phase transition. Hence, in order to prove that the optomechanical model leads to the expected 
phase transition, we have performed simulations of the optomechanical system in the classical limit. Similarly to what 
we did for the DPO model in the first section, this limit is found by assuming that both the optical and mechanical 
modes are in a coherent state at all times. Let us show the procedure explicitly for this case too. Our starting point is 
the original optomechanical master equation but replacing the mechanical dissipator by [x, {p, •}]/2i, where 
p = i{b^ — b) is the mechanical momentum quadrature. For high-Q mechanical oscillators which admit a weak-coupling 
description of their interaction with the environment, this dissipator leads to the same physics as the previous one 
ms], but provides better-looking classical equations. With this change, the evolution of the expectation value of any 
system operator A reads 


d{A) 

dt 


-A 


Af 

dt 


— Hm + ^q(^) + ^qm]) + 7q(([^q: ^]^q) + dq])) + 


7eff 


(53) 


where just as with the DPO model, we are defining the expectation value with respect to the state in the picture 
rotating at the laser frequency, that is, (•) = tr{- p}, with p the state in the rotating frame. Applied to dq, x, and 
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p, and denoting by a{t) and f3{t) the amplitudes of the optical and mechanical coherent states, we find the classical 
evolution equations 

X = Qp, (54a) 

p=-2j^sp-{n-4:gq\af)x, (54b) 

d = £q{t) - (7q - iAq - igqx'^)a, (54c) 

with the classical mechanical position and momentum defined as x = 2Re{/d} and p = 2Im{;d}, respectively. This 
nonlinear system can be efficiently simulated numerically for (practically) any parameter set, and the phonon number 
that it predicts can be evaluated as hmt^oo(^^K^)) = l™t^oo[^^(^) +P^(^)]/4, which again can be approximated by 
n + ^nsin(flqt), with (typically) Sn n. 



